; docformat = 'rst'
;
;+
;
; :Purpose:
;   Write MOZART preprocessor .inp file
;   
; :Inputs:
;   sim_name: name of simulation
;   
; :Requires:
;   <simulation_name>.params file with optional MOZART variables
;
; :Keywords:
;   Flux_output: (Input, Boolean, default=0) set to 1 if fluxes should be output by MOZART
;
; :Outputs:
;   Creates a <simulation_name>.inp file for the MOZART preprocessor in:
;     <cels_input_directory>/<simulation_name>/
;     
; :Example::
;   cels_mozart_preproc, 'SF6', /flux_output
;   
; :History:
; 	Written by: Matt Rigby, MIT, Oct 11, 2011
;
;-
Pro cels_mozart_preproc, sim_name, flux_output=flux_output

  compile_opt idl2
  on_error, 2

  output_name=strlowcase('cels_' + sim_name)
  inp_filename=cels_filestr(/input, sim_name + '/EUM/' + output_name + '.inp')

  species_short=cels_get_parameter(sim_name, 'POLLUTANT_SHORT')
  species_formula=cels_get_parameter(sim_name, 'POLLUTANT_FORMULA')
  Strat_lifetime=cels_get_parameter(sim_name, 'STRAT_LIFETIME')
  OHA=cels_get_parameter(sim_name, 'OHA')
  OHER=cels_get_parameter(sim_name, 'OHER')

  ;Calculate stratospheric lifetime scaling factor
  ;We currently use a stratospheric oxidant field calculated for METHANE
  ;This is simply a scaling for that field
  ;CHANGE THIS FOR YOUR OWN FIELD!!!
  if strat_lifetime ne !null then begin
    print, 'CELS_MOZART_PREPROC: Warning, scaling relative to CH4 stratospheric oxidation!'
    strat_lifetime_scaling=75./strat_lifetime
  endif else begin
    print, 'CELS_MOZART_PREPROC: No stratospheric oxidation'
  endelse
  
  restore, cels_filestr(/Input, sim_name + '/EUM/grid.inp.sav')
  
  emissions_files=file_search(cels_filestr(/input, sim_name + '/emissions/MOZART'), '*')
  
  ;Extract species names from emissions filenames (looks for periods either side of species name)
  species_name=strmid(emissions_files, transpose(strpos(emissions_files, '.')+1))
  species_name=strmid(species_name, 0, transpose(strpos(species_name, '.')))
  
  n_species=n_elements(species_name)
  
  species_formula=replicate(species_formula, n_species)
  
  ;Species name list
  ;Limit to 10 names per line
  n_Species_string=ceil(float(n_species)/10.)
  Species_string=strarr(n_Species_string)
  for si=0, n_species_string-1 do $
    for ti=0, min([n_species - si*10, 10])-1 do Species_string[si]=Strcompress(Species_string[si] + species_name[10*si + ti] + ',')


  ;Begin writing file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  openw, FUN, inp_filename, /get_lun
  
  printf, FUN, 'BEGSIM'
  printf, FUN, ''
  printf, FUN, 'output_file        = ' + output_name + '.doc'
  printf, FUN, 'sim_dat_filename   = ' + output_name + '.dat'
  printf, FUN, ''
  
  printf, FUN, 'Comments'
  printf, FUN, '  Input file for MOZART created by CELS_MOZART_PREPROC'
  printf, FUN, '  ' + sim_name
  printf, FUN, '  ' + systime()
  printf, FUN, 'End Comments'
  printf, FUN, ''
  printf, FUN, 'SPECIES'
  printf, FUN, ''
  printf, FUN, 'Solution'

  ;Write species names and chemical formula
  for si=0, n_species-1 do begin
    printf, FUN, species_name[si], '->', species_formula[si], ','
  endfor
  printf, FUN, 'End Solution'
  printf, FUN, ''
  
  printf, FUN, 'Fixed'
    Fixed_string='M, '
    Fixed_string=Fixed_string + 'OX, '  ;Powerful oxidant used for EUM pollution removal
    if OHA ne !null then Fixed_string=Fixed_string + 'OH, '  ;OH field
    if strat_lifetime ne !null then Fixed_string=Fixed_string + 'SOX, ' ;Stratospheric oxidant
    printf, FUN, Fixed_string
  printf, FUN, 'End Fixed'

  printf, FUN, ''

  printf, FUN, 'End SPECIES'
  printf, FUN, ''

  printf, FUN, 'Solution Classes'
  printf, FUN, ''
  printf, FUN, 'Explicit'
  printf, FUN, species_string
  printf, FUN, 'End Explicit'

  printf, FUN, ''

  printf, FUN, 'End Solution Classes'

  printf, FUN, ''

  printf, FUN, 'CHEMISTRY'
  printf, FUN, ''
  printf, FUN, 'Reactions'

  ;REGULAR REACTIONS
  For si=0, n_species-1 do begin
    if OHA ne !null then begin
      printf, FUN, strcompress(species_name[si] + ' + OH -> products  ;' + $
          string(OHA, format='(E10.2)') + ',' + string(OHER) )
    endif
    if strat_lifetime ne !null then begin
      printf, FUN, strcompress(species_name[si] + ' + SOX -> products  ;' + string(strat_lifetime_scaling, format='(F10.3)'))
    endif
  endfor

  ;Reactions with powerful oxidant
  wh=where(strmatch(species_name, strupcase(strcompress(species_short + 'P*', /remove_all))), count)
  if count gt 0 then begin
    for whi=0, count-1 do begin
      printf, FUN, strcompress(species_name[wh[whi]]  + ' + OX -> products  ; 1.e-21' ) ;THIS REACTION RATE WILL DEPEND ON OXIDANT FIELD USED
    endfor
  endif else begin
    print, 'CELS_MOZART_PREPROC: Warning, no powerful oxidant'
  endelse

  printf, FUN, 'End Reactions'
  printf, FUN, ''
  printf, FUN, 'END CHEMISTRY'
  printf, FUN, ''

  printf, FUN, 'SIMULATION PARAMETERS'
  printf, FUN, ''

  printf, FUN, 'Spatial Dimensions'
  printf, FUN, 'Longitude points = ', string(n_elements(lon))
  printf, FUN, 'Latitude  points = ', string(n_elements(lat))
  printf, FUN, 'Vertical  points = ', string(28)  ;CHANGE FOR NON-NCEP
  printf, FUN, 'End Spatial Dimensions'
  printf, FUN, ''

  printf, FUN, 'Surface Flux'
  ;NAMES HERE
  if keyword_set(flux_output) then for si=0, n_species_string-1 do printf, FUN, species_string[si]
  
  printf, FUN, 'End Surface Flux'
  printf, FUN, ''

  printf, FUN, 'Outputs'
  printf, FUN, ''

  printf, FUN, 'File'
  printf, FUN, 'Transported Species = avrg'
  ;NAMES HERE
  for si=0, n_species_string-1 do printf, FUN, species_string[si]
  printf, FUN, 'End Transported Species'

  printf, FUN, 'End File'
  printf, FUN, ''

  printf, FUN, 'End Outputs'
  printf, FUN, ''

  printf, FUN, 'END SIMULATION PARAMETERS'
  printf, FUN, ''

  printf, FUN, 'ENDSIM'

  close, FUN
  free_lun, FUN
  
End
